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ABSTRACT 

We demonstrate the existence of long-lived, large-amplitude eccentric modes in massive, annular parti- 
cle disks in orbit about a central mass. The lopsided modes we have found precess slowly in the prograde 
direction at a rate which increases with disk mass and decreases with the amplitude of the distortion. 
The lopsidedness generally survives for as long as we run the calculations and may last indefinitely; we 
found no significant decay in over 700 disk-particle orbits in one case. These strongly non-linear modes 
are readily created using a number of different perturbing rules, but we find no evidence for linear insta- 
bilities in unperturbed disks. Our results suggest that Tremaine's eccentric disk model for the nucleus 
of M31 may be viable, but none of our models matches all aspects of the observational data. 

Subject headings: galaxies: nuclei — galaxies: kinematics dynamics — galaxies: Local Group — stars: 
formation 



1. INTRODUCTION 

Massive compact objects are frequently surrounded by 
disks; familiar examples are in accretion disks, protoplan- 
etary disks and planetary ring systems. The disk generally 
has a small mass in comparison with the central object, but 
in some cases it may have sufficient mass to cause a signif- 
icant departure from a Keplerian rotation law. The possi- 
ble small-amplitude, lop-sided (m = 1) modes of such disks 
have been studied by Adams, Ruden & Shu (1989), Shu et 
al. (1990), Tremaine (2000) and others, while Lee & Good- 
man (1999) report a study of possible finite-amplitude, 
one-armed spiral modes. Ogilvie (2001) presents a general 
treatment of non-linear modes in fluid disks. 

Large-amplitude eccentric modes in near-Keplerian par- 
ticle disks have received special attention since Tremaine 
(1995) proposed an eccentric disk model for the nuclear 
structure of the Andromeda galaxy (M31). Sridhar, Syer 
& Touma (1999) and by Salow & Statler (2000) attempt to 
construct self-consistent models, while Bacon et al. (2000) 
have conducted A-body simulations. 

We here report our own A-body simulations of massive, 
eccentric disks in orbit about a central point-like mass. 
We find it surprisingly easy to create long-lived, large- 
amplitude eccentric modes in annular disks having a wide 
range of masses (up to ~ 50% of the central mass), radial 
extents, and eccentricity profiles. Since low-mass disks 
are expected to be cool and thin, we confine the particles 
in our simulations to motion a plane; higher mass disks, 
which are likely to possess larger amounts of random mo- 
tion, may have a non-negligible thickness and require a 3-D 
code. In the final section, we make a preliminary compar- 
ison between our simulations and the data obtained by 
Kormendy & Bender (1999) on the M31 nucleus. 

2. INITIAL CONDITIONS AND METHODS 

Our simulations are of an isolated system with two mass 
components: a razor-thin disk of particles whose motion is 
confined to a plane and a mobile central mass. The central 
object is a slightly softened point mass, M,, having the 
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with a small core radius b. As our initial experiments 
with filled disks did not produce any long-lived modes, 
we here report simulations with annular disks only. Our 
unperturbed disk has the surface density distribution 
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for \R — R \ < a and is zero outside this range. Before 
we introduce eccentricity, the orbits of all particles are 
circular. We adopt M. as our mass unit and i?o as our 
length unit and choose b = 5 x 10~ 3 i?o- We set G — 1, so 
that our time unit is (Rl/GM.) 1 / 2 . The orbital period in 
the middle of the annulus is therefore ~ 2ir. 

In order to facilitate the creation of an initially eccentric 
disk, we begin all our simulations with Md/M, — 10~ 3 , 
and then increase the disk mass to the desired value during 
the first part of the evolution. 

We make a low-mass disk eccentric by adjusting the or- 
bit of each disk particle to be a Kepler ellipse aligned with 
the x-axis and with the central mass at one focus - i.e. 
we neglect the disk contribution to the total potential. 
The particles of the annular disk (eq. 2) are moved and 
given new velocities such that their elliptical orbits inter- 
cept the positive x-axis where their original circular orbit 
would have. Our standard rule (a), sets the distribution 
of eccentricities such that all initial orbits pass through 
the point (x, y) = — {Rq + a, 0) - i.e. the outer edge of the 
disk remains circular while the elliptical inner edge has 
semi-major axis Rq and eccentricity e = a/Ro- We have 
experimented with other rules: (b) making all orbits pass 
through a point outside the original disk so that even the 
outer edge is eccentric, (c) a fixed fraction of the eccentric- 
ity required by these rules, (d) a disk that is eccentric only 
in the inner parts and axisymmetric farther out, (e) a con- 
stant eccentricity throughout the disk, and (f) a reversed 
trend with radius so that the inner edge is circular. 
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Fig. 1. — A short sequence of snapshots from our canonical 
model (only 1 particle in 20 is plotted in each frame). The 
position of the central mass is marked by the plus symbol. 



We allow the low-mass disk to evolve for a few orbital 
periods, although little happens during this time. We then 
raise the mass of each particle so that the disk mass rises at 
a constant rate. Generally, this rate is dM^/dt = 1CP 3 M., 
so that it takes ~ 16 orbital periods for the disk mass to 
rise by 0.1M.; other rates were also tried. Once the desired 
mass is reached, we reposition our coordinate origin at the 
center of mass and subtract any net momentum from the 
system. The subsequent evolution is unrestrained. 

We evaluate gravitational forces between the disk par- 
ticles using a grid method for efficiency (Sellwood 1997). 
A 2-D Cartesian grid (e.g. Hohl & Hockney 1969) is ideal 
for this case because the disk is thin and not centrally 
concentrated, and the grid has no preferred center. We 
determine the forces between the central mass and disk 
particles directly and compute the motion of the central 
mass in response to the total force from the disk. We sub- 
divide the time steps of disk particles as they approach the 
central object. The code conserves energy and both linear 
and angular momentum to a high degree of precision. We 
employ 100 000 disk particles, a 256 2 grid, and a time step 
of 0.01 or 0.005. We have verified that our results are not 
significantly different when the grid size, particle number 
or time step is varied by factors of two or more. 

3. RESULTS 

Our canonical model has a final disk mass of 0.1M., 
a = Rq/2 and the initial eccentricity prescription (a) above 
(i.e. the ratio of apocenter to pericenter radius is three for 
the innermost orbit). As the disk mass was increased, the 
initially neat elliptical stream pattern developed transient 
spiral patterns, heated somewhat and spread radially. As 
soon as the disk mass stopped increasing (at t = 200), 
the model settled to a lop-sided disk with no significant 
spirality, as shown in Figure |l|. 
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Fig. 2. — Variation of pattern speed with disk to central mass 
ratio. The best linear fit is shown by the line. 



The major axis of the eccentric disk precessed at an 
almost steady rate with no obvious decay. To quantify 
this behavior, we form the Fourier coefficients 
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where (Rj, 4>j) are polar coordinates of the j-th particle at 
time t. The amplitude of the m = 1 distortion decreased 
slightly from its initial value as the disk relaxed while the 
mass was increased. Subsequently, the amplitudes of the 
m = 1 coefficients manifest a mild modulation, suggesting 
beats between two waves with a gradually increasing beat 
period. The frequencies (pattern speeds) of the stronger 
and weaker waves are Vt p ~ 0.045 and f2 p ~ 0.008 in our 
units, in the first part of the evolution, but the stronger 
wave later slows to Cl p ~ 0.040. The amplitude of neither 
wave decays significantly during the last ~ 4500 dynamical 
times - these disturbances appear to be indefinitely-lived, 
large-amplitude, lop-sided modes. 

We were able to find long-lived modes in most of the 
annular-disk models we tried. Models with a = Q.5Rq 
and a higher disk mass behaved reasonably well for Md 5s 
0.3M.. The large random motions in disks of yet higher 
mass cause small numbers of particles to become unbound; 
the surviving disk appears to remain lop-sided, but the 
drift in the center of mass of the bound particles com- 
plicates the analysis. The pattern speed of the domi- 
nant mode, which is always the fastest, increases in an 

f proximately linear manner with the disk mass (Figure 
, behavior that was predicted by Tremaine (2000) for 
small-amplitude waves. A linear fit yields fi p ~ (0.393 ± 
0.008)Md/M». We stress that our results are for waves 
of considerable amplitude and we report below that less 
strongly perturbed disks have higher pattern speeds, in 
agreement with Salow & Statler (2000). 

Reducing the initial eccentricity, but holding a = Rq/2 
and Md/M, = 0.1, increased the pattern speed of the re- 
sulting mode but had little effect on its survivability, at 
least until the distortion was very mild. A lopsided mode 
stood out clearly in a model in which the initial eccentricty 
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of all orbits was reduced to 20% of that required for rule 
(a), but a long-lived mode was barely detectable when the 
initial distortion was reduced to 10%. 

In fact, we were able to generate long-lived eccentric 
disks with quite remarkable ease. Lopsidedness survived, 
but usually with different pattern speeds, in all disks hav- 
ing strong initial eccentricities, for almost all types of ini- 
tial perturbation. We also varied the radial extent of the 
disk in the range 1/2 < a/Ro < 7/9; low-mass disks settled 
to a long-lived mode in every case. Cases in which the lop- 
sidedness weakened considerably were for more massive 
disks (<; 0.3M.) with larger radial extents (a/Ro ^ 2/3), 
whether the initial perturbation included all or just part 
of the disk. 

Statler (1999) defines orbital eccentricity e = (x_ + 
x + )/(.t_ — x + ), where x+ and a;_ are respectively the in- 
tercepts of an aligned periodic orbit on the positive and 
negative x-axis. He argues that a self-consistent mode re- 
quires e to change sign near the densest part. We use the 
positions and velocities of the particles at one instant to 
define Kepler ellipses for each. As they are, in general, 
misaligned with the principal axis of the mode, we define 
e from the intercepts with the ellipse major axis. The disk 
appears to have two populations of particles, a strongly 
eccentric group which all have e ^ 0.7 and are generally 
aligned with the mode, and a "hot" population of particles 
with —0.4 ^ e ^ 0.4 whose orbits are aligned randomly 
with respect to the major axis of the mode. This structure 
appears to be generic, although the relative proportions of 
the two populations in different models vary. Thus, we 
find no evidence of a coherent change in orientation in the 
dense part of the disk. 

We were unable to detect any growing modes in a disk 
with no initial distortion; transient spiral patterns as the 
disk mass increased heated the particle distribution to 
Q r*> 3 in the densest part of the disk, but it then re- 
mained closely axisymmetric. Tremaine (2000) predicts 
that small-amplitude modes should be stable, and our re- 
sults also suggest that no lopsided instability is present. 

4. COMPARISON WITH M31 

The brightest point of M31 is displaced slightly from the 
center of the bulge isophotes (Light et al. 1974) and the nu- 
cleus is now known to be double (Lauer et al. 1993), with 
two brightness peaks separated by ~ 0.5". The brighter 
peak is designated PI, while the fainter (P2) lies within 
0.05" of the center of the outer nuclear isophotes and 
the bulge photocenter. Ground-based integral-field spec- 
troscopy by Bacon et al. (1994) revealed that P2 is close 
to the rotational center of the system, and that the dis- 
persion peak coincides with neither brightness peak, and 
is displaced from P2 in a direction almost exactly oppo- 
site to PI. Kormendy & Bender (1999, hereafter KB99) 
confirm these features from spectroscopy with a single-slit 
with higher-resolution. Bacon et al. (2000) report yet more 
detailed observations which confirm and extend previous 
work and which we discuss in a later paper. 

In order to account for the observations, Tremaine 
(1995) proposed a model in which an eccentric disk of 
stars orbits a black hole (BH). The fainter peak, P2, is 
presumed to be the location of the BH at the center of the 
large spheroidal bulge, while the brighter peak is caused 
by the lingering of stars at apocenter in the eccentric disk, 
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Fig. 3. — Comparison of the KB99 data with some of our projected 
models. The observed kinematics along a cut passing through PI 
and P2 are shown by crosses with error bars in panels (b) and (c) 
only; PI is to the right of the origin and KB99 estimate that the BH 
(P2) is 0.080 ± 0.011" to the left of the rotation center. Our best 
model is shown by the full-drawn lines and is described in the text. 
The disk masses in the other models are M^/M, = 0.001 (dotted), 
0.05 (dashed), 0.1 (dot-dashed) and 0.3 (triple-dot-dashed). The 
radial extents of the disks and the imposed perturbations were all 
different and will be described in a later paper. 

viewed almost edge-on. His model, which was designed 
to fit the photometric asymmetries only, placed stars in 
the disk on nested confocal elliptical orbits around the 
BH which dominates the total mass of the nucleus. KB99 
claim general agreement between their data and the kine- 
matic predictions of Tremaine's model in both the shape 
and asymmetry of the circular speed profile about P2. 

We here attempt a quantitative comparison between our 
self-consistent dynamical models and the KB99 data. The 
symmetry in the bulge rotation curve well outside 3" al- 
lows KB99 to determine the systemic velocity of M31 to a 
precision of better than 1 km s _1 . They adopt a slightly 
different photometric bulge model from that proposed by 
Tremaine, but in either case the circular speeds arising 
from the central attraction of the bulge for M/Ly = 4 
are negligible in comparison with those from the BH and 
nuclear disk in the inner 2". KB99 subtract a scaled spec- 
trum of the outer bulge from their spectra in the nuclear 
region to yield "bulge subtracted" estimates of the nuclear 
disk kinematics. These data indicate streaming speeds ris- 
ing to ~ 230 km s _1 on the PI side, but to only ~ 180 km 
s _1 on the other side [see Figure 0(c)]. The asymmetry 
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in velocities persists out to <; 3" from the rotation center, 
which they find is very close to P2 (0.080 ±0.011" towards 
PI). 

Although Tremaine suggests that the eccentric nuclear 
disk is coplanar with the outer disk of M31 (i.e. inclined 
at 77° from face-on), we "view" our razor-thin disk models 
edge-on in order to avoid having to thicken our model in 
some artificial manner. We rotate our models so that the 
azimuthal angle in the disk-plane between disk major axis 
and the plane of the sky is 21°, as deduced by Tremaine. 
We then compute the line of sight velocity and dispersion 
from the projected particle speeds using a kernel (Silver- 
man 1986) with a fixed width chosen to approximate the 
spatial resolution of the KB99 data (ct* = 0.27"). We 
fit each of our models to the observed streaming velocity 
profile to ±1" only, fitting the projected scales in both dis- 
tance and velocity, together with a possible position offset, 
to minimize \ 2 ■ Note that we do not allow the systemic 
velocity to be a free parameter in the fits. 

A characteristic signature of the lopsided disk is the 
asymmetry in the streaming velocity on the two sides of 
the nucleus. We have therefore focused on trying to fit the 
streaming velocity pattern. Our most successful model 
is shown by the full-drawn curve in Figure |(c). This 
model has M d /M. = 0.2, a/R = 2/3 and both the in- 
ner and outer edges were initially elliptic with all orbits 
passing through the point (— 2i?o,0). The evolved model 
(full-drawn curve) reproduces the ~ 50 km s _1 difference 
between the peak velocities observed on the two sides tol- 
erably well, and the streaming velocities in the model de- 
crease from the central peaks roughly as do the data. The 
full-drawn curves in panels (a) and (b) of Figure || show, 
respectively, the projected mass density and velocity dis- 
persion of the same model. The mass density peaks at 
about 0.6" from the center, slightly farther than the ob- 
served position of PI (0.5"); the peak value is ~ 0.6 mag 
above that at P2 which is nicely consistent with the cut 
through the HST photometry shown by KB99. Also the 
position of the BH, marked by the vertical line, shows a 
small shift from the rotation center towards P2 of about 
the amount estimated by KB99. The velocity dispersion 
profile does not fit the KB99 data at all well, however; the 
model is too cool everywhere and the position of the peak, 
while on the correct side, is too far from the center. 

We have tried a large number of other models, but none 
fits the data as well. Some of the better ones are shown by 
the broken curves in Figure |[ A disk of virtually massless 
particles (dotted curves) with a broad radial extent has too 
great an asymmetry and the projected velocity dispersion 
is the lowest of all. We have found that much heavier disks 
(e.g. Md/M, = 0.3, the triple-dot-dashed curve) tend not 
to show strong enough asymmetry in the streaming speed 
and the velocity dispersion is actually lower in this case 
than in our preferred model. 

The scaling of our preferred model to the KB99 data 
require i?o = 0.523" = 1.95 pc (for an adopted distance of 
770 kpc) and a velocity unit = 301.6 km s^ 1 . With this 
scaling, one dynamical time is approximately 6 300 years 
and the implied BH mass is M. ~ 4.0 x 10 7 M Q . 

5. CONCLUSIONS 

We have shown that lopsided modes in an almost- 
Keplerian annular disk can survive for at least 700 dynami- 



cal times, and probably indefinitely. The modes, which are 
easily excited with a range of amplitudes in disks of dif- 
fering radial extents, survive with significant amplitude to 
disk masses at least as large as 30% of the central mass. If 
the modes are also present in collisional, or gaseous disks, 
they could be of relevance to the formation of massive 
planets or small companion stars from a disk surrounding 
a primary, or for the persistence of asymmetric accretion 
disks in other contexts. 

The existence of these long-lived modes offers consider- 
able support to Tremaine's (1995) eccentric disk model for 
the double nucleus of M31, but we have been unable to cre- 
ate a model that matches all aspects of the observations. 
One of our models reproduces some significant kinematic 
and photometric features of the M31 nucleus quite well, 
including the off-center brightness peak, the asymmetry 
in the streaming velocities, and approximately correct lo- 
cation of the velocity zero between PI and P2. The veloc- 
ity dispersion peak is also off-center, but not in quite the 
right place and the model is altogether too cool. Tremaine 
has suggested to us that two-body relaxation is important 
in the nucleus, a process largely absent from our models. 
It remains to be seen whether a model with the consid- 
erably greater random motion needed to fit the observed 
dispersion profile could continue to fit the streaming and 
photometric profiles. 

We thank Scott Tremaine for his interest and advice 
and the referee for a thoughtful report. This work was 
supported by NASA LTSA grant NAG 5-6037. 
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